# This file is part of Sonedyan.
#
# Sonedyan is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public
# License as published by the Free Software Foundation;
# either version 3 of the License, or (at your option) any
# later version.
#
# Sonedyan is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied
# warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
# PURPOSE.  See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public.
# If not, see <http://www.gnu.org/licenses/>.
#
# Copyright (C) 2009-2012 Jimmy Dubuisson <jimmy.dubuisson@gmail.com>

source("misc.R")

# http://snipplr.com/view/9914/
#def assortativity(graph, degrees=None):
#      	if degrees is None: degrees = graph.degree()
#      	degrees_sq = [deg**2 for deg in degrees]
#
#      	m = float(graph.ecount())
#      	num1, num2, den1 = 0, 0, 0
# 
#    	for source, target in graph.get_edgelist():
#      		num1 += degrees[source] * degrees[target]
#      		num2 += degrees[source] + degrees[target]
#      		den1 += degrees_sq[source] + degrees_sq[target]
#
#      	num1 /= m
#      	den1 /= 2*m
#      	num2 = (num2 / (2*m)) ** 2
#
#      	return (num1 - num2) / (den1 - num2)
assortativity <- function(g)
{
	degrees <- degree(g, V(g), mode = "out")
	degrees_sq <- degrees ** 2
	
	m <- length(E(g))
	num1 <- 0
	num2 <- 0
	den1 <- 0
	edges <- get.edgelist(g, names = FALSE)
	
	for (i in c(1:m))
	{
		source <- as.numeric(edges[,1][i]) + 1
		target <- as.numeric(edges[,2][i]) + 1
		
		num1 <- num1 + (degrees[source] * degrees[target])
		num2 <- num2 + (degrees[source] + degrees[target])
		den1 <- den1 + (degrees_sq[source] + degrees_sq[target])
	}
	
	num1 <- num1 / m
	den1 <- den1 / (2 * m)
	num2 <- (num2 / (2 * m)) ** 2
	
	(num1 - num2) / (den1 - num2)
}

# compute graph efficiency
efficiency <- function(g)
{
	n <- length(V(g))
	mat <- shortest.paths(g, V(g), mode = "out")
	
	s <- 0
	
	for (i in c(1:(n - 1)))
	{
		for (j in c((i + 1):n))
		{
			s <- s + (1 / mat[i, j])
		}
	}
	
	s / (n * (n - 1))
}

# get the avg shortest path & size of the greatest 
# connected component of the specified directed graph
#
# g: graph to be measured
# ml: list of metrics to be returned
#
# cpl: characteristic path length
# msccs: main SCC size
# mwccs: main WCC size
# ddiam: directed diameter
# udiam: undirected diameter
# eff: efficiency
# assort: assortativity
# gtrans: global transitivity (global CC)
# altrans: avg of local transitivities (avg of local CCs)
# dens: density
# nv: # vertices
# ne: # edges
# avgpr: avg pagerank
# avgdeg: avg degree
# avgindeg
# avgoutdeg
# maxdeg
# maxindeg
# maxoutdeg
# mindeg
# minindeg
# minoutdeg 
get.graph.stats <- function(g, ml)
{
	results <- c()

	for (m in ml)
	{
		if (m == "cpl")
			results["cpl"] <- average.path.length(g)
		else if (m == "msccs")
			results["msccs"] <- get.main.cc.size(g)
		else if (m == "mwccs")
			results["mwccs"] <- get.main.cc.size(g, mode = "weak")
		else if (m == "ddiam")
			results["ddiam"] <- diameter(g, directed = TRUE)
		else if (m == "udiam")
			results["udiam"] <- diameter(g, directed = FALSE)
		else if (m == "eff")
			results["eff"] <- efficiency(g)
		else if (m == "assort")
			results["assort"] <- assortativity(g)
		else if (m == "gtrans")
			results["gtrans"] <- transitivity(g, type = "global")
		else if (m == "altrans")
			results["altrans"] <- mean(transitivity(g, vids = V(g), type = "local"))
		else if (m == "dens")
			results["dens"] <- graph.density(g, loops = FALSE)
		else if (m == "nv")
			results["nv"] <- length(V(g))
		else if (m == "ne")
			results["ne"] <- length(E(g))
		else if (m == "avgpr")
			results["avgpr"] <- mean(page.rank(g, vids = V(g))$vector)
		else if (m == "avgdeg")
			results["avgdeg"] <- mean(degree(g, mode = "total"))
		else if (m == "avgindeg")
			results["avgindeg"] <- mean(degree(g, mode = "in"))
		else if (m == "avgoutdeg")
			results["avgoutdeg"] <- mean(degree(g, mode = "out"))
		else if (m == "maxdeg")
			results["maxdeg"] <- max(degree(g, mode = "total"))
		else if (m == "maxindeg")
			results["maxindeg"] <- max(degree(g, mode = "in"))
		else if (m == "maxoutdeg")
			results["maxoutdeg"] <- max(degree(g, mode = "out"))
		else if (m == "mindeg")
			results["mindeg"] <- min(degree(g, mode = "total"))
		else if (m == "minindeg")
			results["minindeg"] <- min(degree(g, mode = "in"))
		else if (m == "minoutdeg")
			results["minoutdeg"] <- min(degree(g, mode = "out"))
	}
	
	results
}

##########
## display functions
##########

# display max degree vertices
#
# g: graph
# dd: graph degrees
# m: mode ("in", "out", "tot")
# r: rank 
display.max.degree.vertices <- function(g, dd, m, r)
{
	md <- which(dd != 0) - 1
	mdd <- md[length(md) - (r - 1)]
	mdids <- V(g)[(which(degree(g, mode = m) == mdd)) - 1]
	print(paste("Max (rank ", r, ") ", m, "-degree (", mdd, ") vertices: ", V(g)[mdids]$name, sep = ""))
}

# print global stats of the specified graph
print.graph.global.stats <- function(g, name)
{
	print(paste("Printing stats for '", name, "'"))
	print(paste("# vertices: ", length(V(g))))
	print(paste("# edges: ", length(E(g))))
	print(paste("Density: ", graph.density(g, loops = FALSE)))
	print(paste("Diameter (directed): ", diameter(g, directed = TRUE)))
	print(paste("Diameter (undirected): ", diameter(g, directed = FALSE)))
	print(paste("Transitivity (CC): ", transitivity(g)))
	print(paste("Characteristic path length: ", average.path.length(g)))
	print(paste("Efficiency: ", efficiency(g)))
	print(paste("Assortativity: ", assortativity(g)))
	
	din <- degree(g, mode = "in")
	dout <- degree(g, mode = "out")
	dtot <- degree(g, mode = "total")
	
	print(paste("Avg / Max / Min / Var in-degree: ", mean(din), "/", max(din), "/", min(din), "/", var(din)))
	print(paste("Avg / Max / Min / Var out-degree: ", mean(dout), "/", max(dout), "/", min(dout), "/", var(dout)))
	print(paste("Avg / Max / Min / Var degree: ", mean(dtot), "/", max(dtot), "/", min(dtot), "/", var(dtot)))
	
	ddin <- degree.distribution(g, mode = "in")
	ddout <- degree.distribution(g, mode = "out")
	ddtot <- degree.distribution(g, mode = "total")
	
	#print(paste("Max in-degree vertices: ", V(g)[(which(degree(g, mode = "in") == max(din)) - 1)]$id))
	#print(paste("Max out-degree vertices: ", V(g)[(which(degree(g, mode = "out") == max(dout)) - 1)]$id))
	#print(paste("Max degree vertices: ", V(g)[(which(degree(g, mode = "total") == max(dtot)) - 1)]$id))
	
	mdin <- which(ddin != 0) - 1
	mdout <- which(ddout != 0) - 1
	mdtot <- which(ddtot != 0) - 1
	
	print("----- Max in-degree vertices -----")
	
	for (i in c(1:10))
		display.max.degree.vertices(g, ddin, "in", i)
	
	print("----- Max out-degree vertices -----")
	
	for (i in c(1:5))
		display.max.degree.vertices(g, ddout, "out", i)
	
	print("----- Max degree vertices -----")
	
	for (i in c(1:10))
		display.max.degree.vertices(g, ddtot, "total", i)
}

# print all stats of the specified graph
print.graph.all.stats <- function(g, name)
{
	print(paste("Printing stats for '", name, "'"))
	print(paste("# vertices: ", length(V(g))))
	print(paste("# edges: ", length(E(g))))
	print(paste("Density: ", graph.density(g, loops = FALSE)))
	print(paste("Diameter (directed): ", diameter(g, directed = TRUE)))
	print(paste("Diameter (undirected): ", diameter(g, directed = FALSE)))
	print(paste("Transitivity (CC): ", transitivity(g)))
	print(paste("Characteristic path length: ", average.path.length(g)))
	print(paste("Efficiency: ", efficiency(g)))
	print(paste("Assortativity: ", assortativity(g)))
	
	# get graph clusters
	wccs = clusters(g, mode="weak")
	sccs = clusters(g, mode="strong")
	
	print(paste("# wccs: ", wccs$no))
	print(paste("# sccs: ", sccs$no))
	
	print(paste("max wcc: ", max(wccs$csize)))
	print(paste("max scc: ", max(sccs$csize)))
}

#
# nv: # vertices
# ne: # edges
# steps: # of steps
# cl: length of cycles
get.er.metrics <- function(nv, ne, steps, cl)
{
	# ae: avg efficiency
	# at: avg transitivity
	# aa: avg assortativity
	gmetrics <- hash()
	gmetrics["ae"] = 0
	gmetrics["at"] = 0
	gmetrics["aa"] = 0
	
	# avcov: avg vertices coverage
	# aecov: avg edges coverage
	# anc: avg # of cycles
	cmetrics <- hash()
	cmetrics["avcov"] = 0
	cmetrics["aecov"] = 0
	cmetrics["anc"] = 0
	
	for (i in c(1:steps))
	{
		#er <- erdos.renyi.game(nv, d, type = "gnp", directed = TRUE)
		er <- erdos.renyi.game(nv, ne, type = "gnm", directed = TRUE)
		V(er)$id <- c(1:length(V(er)))
		
		cs <- get.cycles(er, cl)
		
		vcov <- length(cs$v) / nv
		cmetrics["avcov"] <- cmetrics[["avcov"]] + vcov
		ecov <- length(cs$e) / ne
		cmetrics["aecov"] <- cmetrics[["aecov"]] + ecov
		nc <- length(cs$c)
		cmetrics["anc"] <- cmetrics[["anc"]] + nc
		
		eff <- efficiency(er)
		gmetrics["ae"] <- gmetrics[["ae"]] + eff
		
		trans <- transitivity(er)
		gmetrics["at"] <- gmetrics[["at"]] + trans
		
		ass <- assortativity(er)
		gmetrics["aa"] <- gmetrics[["aa"]] + ass
		
		#loginfo(paste("Vertices coverage: ", vcov))
		#loginfo(paste("Edges coverage: ", ecov))
		#loginfo(paste("# cycles: ", nc))
		
		#loginfo(paste("graph efficiency: ", eff))
		#loginfo(paste("graph transitivity: ", trans))
		#loginfo(paste("graph assortativity: ", ass))
	}
	
	# normalizing
	cmetrics["avcov"] <- cmetrics[["avcov"]] / steps
	cmetrics["aecov"] <- cmetrics[["aecov"]] / steps
	cmetrics["anc"] <- cmetrics[["anc"]] / steps
	
	gmetrics["ae"] <- gmetrics[["ae"]] / steps
	gmetrics["at"] <- gmetrics[["at"]] / steps
	gmetrics["aa"] <- gmetrics[["aa"]] / steps
	
	list(g = gmetrics, c = cmetrics)
}

# printing ER metrics
print.er.metrics <- function(metrics)
{
	cmetrics <- metrics$c
	gmetrics <- metrics$g
	
	print(paste("Avg vertices coverage: ", cmetrics[["avcov"]]))
	print(paste("Avg edges coverage: ", cmetrics[["aecov"]]))
	print(paste("Avg # cycles: ", cmetrics[["anc"]]))
	
	print(paste("Avg graph efficiency: ", gmetrics[["ae"]]))
	print(paste("Avg graph transitivity: ", gmetrics[["at"]]))
	print(paste("Avg graph assortativity: ", gmetrics[["aa"]]))
}

# get colink coefficient
#
# g: graph
# cog: colink graph
# i: vertex id
get.colink.coefficient <- function(g, cog, i)
{
	di <- degree(g, i, mode = "out")
	
	vid <- V(g)[i]$id
	i2 <- which(V(cog)$id == vid) - 1
	
	# if it is numeric(0)
	if (length(i2) == 0)
		return(0)	
	else
		return(degree(cog, i2, mode = "total")/di)
}

# get local clustering coefficient of specified node 
get.clustering.coefficient <- function(g, i)
{
	trans <- transitivity(g, i, type = "local")
	
	if (is.nan(trans))
		return(0)
	else
		return(trans)
}

# get page ranks of all graph nodes
get.page.ranks <- function(g)
{
	page.rank(g, V(g))$vector
}

# get robustness
# NB: measures the proportion of the graph that remains strongly connected under node attacks (averaged over all nodes)
get.robustness <- function(g, mode = "strong", threshold = 1)
{
	l <- length(V(g))
	# avg robustness
	a <- 0
	# list of 'weak' vertices
	wp <- c()
	
	for (i in V(g))
	{
		ng <- delete.vertices(g, i)
		nvids <- get.main.cc.vertex.indexes(ng, mode)
		lr <- length(nvids) / (l - 1)
		
		if (lr < threshold)
		{
			wp <- append(wp, i)
		}
		
		a <- a + lr
	}
	
	return(list(r = a/l, w = wp))
}

# get network entropy as defined in 
# "Robustness and network evolution - An entropic principle"
get.network.entropy <- function(adj, normalize = FALSE)
{
	trans <- get.pev.transition.matrix(adj)
	stat <- get.stationary.distribution(trans)
	er <- get.entropy.rate(trans, stat)
	
	if (normalize)
	{
		return(er / log(length(stat)))
	}
	else
	{
		return(er)
	}
}

# get entropy of degree distribution
get.degree.entropy <- function(g, normalize = FALSE)
{
	distr <- degree.distribution(g)
	
	if (normalize)
	{
		return(get.entropy.2(distr) / log(length(V(g))))
	}
	else
	{
		return(get.entropy.2(distr))
	}
}